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SUMMARY 

The  next  generation  of  commercial  aircraft  will  include 
turbofnn  engines  with  performance  levels  significantly  better 
than  those  in  the  current  fleet.  Control  of  particulate  and 
gaseous  emissions  will  also  be  an  integral  part  of  the  engine 
design  criteria.  These  performance  and  emission  requirements 
present  a  technical  challenge  for  the  combustor:  control  of  the 
fuel  and  air  mixing  and  control  of  the  local  stoichiometry  will 
have  to  be  maintained  much  more  rigorously  than  with  com¬ 
bustors  in  current  production.  A  better  understanding  of  the 
flow  physics  of  liquid  fuel  spray  combustion  is  necessary. 
This  paper  describes  recent  experiments  on  spray  combustion 
where  detailed  measurements  of  the  spray  characteristics  were 
made,  including  local  drop-size  distributions  and  velocities. 
Also,  an  advanced  combustor  CFD  code  has  been  under  deve¬ 
lopment  and  predictions  from  this  code  are  compared  with 
experimental  results.  Studies  such  as  these  will  provide 
information  to  the  advanced  combustor  designer  on  fuel  spray 
quality  and  mixing  effectiveness.  Validation  of  new  fast, 
robust,  and  efficient  CFD  codes  will  also  enable  Ihe  combus¬ 
tor  designer  to  use  them  as  valuable  additional  design  tools  for 
optimization  of  combustor  concepts  for  the  next  generation  of 
aircraft  engines. 

A.  INTRODUCTION 

Aircraft  engines  being  envisioned  now  for  the  next-century 
aircraft  will  have  requirements  that  present  formidable  techni¬ 
cal  challenges  to  the  combustor  designer.  In  the  subsonic 
commercial  transport  arena,  demand  for  low  operating  cost 
translates  into  reduced  fuel  consumption  and  improved  dura¬ 
bility  and  reliability.  Higher  operating  pressures  of  the 
combustor  are  forecast,  with  higher  resulting  fuel  system 
turndown  ratios.  Coupled  with  these  performance  requirements 
is  the  demand  world  wide  for  control  of  pollutant  emissions 
from  aircraft  engines,  especially  oxides  of  nitrogen,  in  the 
arena  of  supersonic  commercial  transports,  forecasts  are 


predicting  fleets  of  hundreds  of  transports  operating  at 
Mach  2.0  to  2.5,  with  airfares  only  slightly  higher  than  today's 
long-range  subsonic  fares.  Economical  fuel  consumption  is  a 
requirement.  But  the  largest  technical  challenge  here  is  low 
levels  of  oxides  of  nitrogen  emissions  during  engine  cruise 
such  that  there  would  be  no  adverse  impact  to  the  earth's 
environment,  specifically  the  atmosphere's  ozone  layer. 

These  technical  and  environmental  challenges  represent  design 
requirements  for  the  combustor  that  are  outside  of  engine 
companies'  experience.  Empirically  based  design  methods  are 
insufficient  by  themselves.  To  augment  this  design  system,  the 
companies  are  increasingly  turning  to  computational  fluid 
dynamics  (CFD)  computer  codes.  Severe  limitations  with  cur¬ 
rently  available  codes  are  excessive  time  requirements  to  run 
a  CFD  code  to  analyze  a  complex  combustor  design:  a  penalty 
both  in  turn-around  time  for  design  answers  as  well  as  cost  for 
the  calculations.  We  at  NASA  are  attempting  to  improve  this 
current  situation  in  the  industry  by  developing  a  fast,  robust, 
efficient  computer  code  for  internal  chemical  reacting  flows. 
The  objective  is  to  produce  a  CFD  code  that  can  be  used  as 
a  more  powerful  design  tool  in  the  industry  to  analyze  com¬ 
plex  combustor  designs  in  significantly  shorter  lum-around 
time  than  current  computer  codes  can  achieve. 

An  integral  pari  of  Ihe  development  of  a  new  CFD  code  is  the 
validation  of  this  code  with  experiments  that  represent  the 
complex  features  of  the  flows  which  need  to  be  analyzed. 
Therefoie.  at  NASA  we  are  also  conducting  experiments  on 
liquid  fuel  spraj  combusting  flows  with  increasingly  complex 
features.  This  will  provide  some  of  the  required  data  for 
validation  of  our  CFD  code  and  other's  codes.  Detailed  mea¬ 
surements  of  both  liquid  sprays  and  gas  characteristics  are 
being  obtained  under  both  nonburning  and  burning  conditions. 

This  paper  will  highlight  recent  results  that  NASA  has 
obtained  with  its  spray  combustion  experiment  and  describe 
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the  formulation  and  physical  modeling  of  its  new  spray  com¬ 
bustion  CFD  code  (ALLSPD).  The  application  of  this  code  to 
combustion  problems  will  be  illustrated  by  several  examples. 

B.  SPRAY  COMBUSTION  EXPERIMENTS 

Combusting  sprays  are  very  important  for  gas  turbine  engine 
applications.  The  investigation  of  combusting  spiays  should 
lead  to  a  better  understanding  of  die  physics  involved  in  this 
complicated  process.  Important  processes  involved  in  combust¬ 
ing  sprays  are  the  interactions  between  the  droplets  and  the 
gas  phase,  the  vaporization  of  the  droplets,  and  chemical 
reactions  with  heat  release.  These  physical  processes  are 
coupled  and  can  only  be  completely  described  using  numerical 
modeling.  As  part  of  an  effort  to  improve  the  numerical 
modeling  of  gas  turbine  combustors,  an  experimental  study 
has  been  undertaken  to  obtain  a  data  set  for  a  relatively  simple 
liquid-fueled  combustor  lliai  can  be  used  for  comparison  widi 
numerical  models  (Ref.  I). 

Because  of  their  practical  applications,  swirling  flows  with 
combustion  have  been  studied  by  a  large  number  of  investi¬ 
gators.  Earlier  reviews  of  swirling  flows  both  with  and  without 
combustion  present  some  general  trends  (Refs.  2  to  4).  These 
papers  predate  die  development  of  nonintrusive.  laser-based 
diagnostics;  consequently  all  of  the  results  described  were 
obtained  using  intrusive  instrumentation  and  detailed  structure 
measurements  for  these  types  of  flows  were  not  possible.  With 
the  advent  of  newer  instnimentation  techniques,  namely  laser 
Doppler  anemometry.  additional  details  of  die  structure  of 
these  types  of  flows  began  to  emerge.  Laser  Doppler  aneino- 
metry  velocity  measurements  in  spray  flames  (Refs.  5  and  6) 
reveal  some  of  the  flowfield  structure  of  swirling  flames.  The 
development  of  the  phase/Doppler  particle  analyzer  (Ref.  7). 
enabled  the  simultaneous  measurement  of  droplet  size  and 
velocity.  This  instrument  has  been  used  by  a  number  of  inves¬ 
tigators  for  measurements  in  spray  flames  in  a  variety  of 
configurations  (Refs.  8  to  12).  This  instrument  has  the 
capability  to  measure  velocities  of  both  the  gas  and  droplet 
phases  in  a  combusting  spray. 

B.l  EXPERIMENT 

The  combustor  utilized  in  the  present  experiment  is  illustrated 
in  Fig.  I.  It  consists  of  a  center  mounted  air-assist  fuel  nozzle, 
Parker  Hannifin  research  simplex  air-assist  atomizer,  sur¬ 
rounded  by  a  coflowing  air  stream.  The  nozzle  orifice  diame¬ 
ter  was  4.8  mm.  Both  the  air  assist  and  die  coflow  air  streams 
had  swirl  imparted  to  them  using  45  degree  swirlers.  The 
swirlers  were  constructed  by  machining  45  degree  slots  into 
rings.  Both  streams  were  swirled  in  the  same  direction  for  the 
present  study.  The  combustion  air  was  not  preheated  and 
entered  the  combustor  at  297  K.  The  top  of  the  air-assist 
nozzle  was  water  cooled  to  prevent  overheating  of  an  O-ring 
in  die  nozzle  assembly.  The  temperatures  of  the  fuel,  atomiz¬ 
ing  air  and  coflow  air  streams  were  measured  using  Chroma! 
Alumel  thermocouples.  Row  rates  of  the  air  streams  were 
measured  using  calibrated  orifices  and  the  fuel  flow  rate  was 
measured  using  a  mass  flowmeter.  All  results  reported  in  the 
present  study  are  reported  for  a  coflow  air  flow  rate  of 
13.88  g/s,  an  air-assist  flow  rate  of  0.96  g/s.  and  a  fuel  flow 
rate  of  0J8  g/s  The  fuel  used  was  heptane.  The  coflow 


stream  entered  the  combustor  in  diree  radial  locations,  passed 
through  a  honeycomb  flow  straightener,  and  the  swirlers 
before  exiting  the  combustor.  The  swirler  was  located  140  mm 
upstream  of  the  combustor  exit.  The  flow  from  the  combustor 
discharged  into  ambient,  stagnant  surroundings. 

The  combustor  was  mounted  vertically  within  a  large 
(1.8  m  sq  by  2.4  in  high)  enclosure.  The  entire  enclosure  was 
mounted  on  two  sets  of  linear  bearings  and  was  traversed 
using  stepper  motors  to  provide  motion  in  two  directions.  The 
combustor  assembly  itself  could  be  traversed  in  the  vertical 
direction  using  a  diird  stepper  motor  to  allow  measurements 
at  all  locations  in  the  flowfield.  This  arrangement  allowed 
rigid  mounting  of  all  optical  components. 

The  phns e/Doppler  particle  analyzer  was  used  for  all  mea¬ 
surements  reported  in  this  study.  A  schematic  of  the  two- 
component  instrument  is  shown  in  Fig.  2.  The  beam  from  a 
6  W  Argon-Ion  laser  is  split  into  488.0  and  514.5  mn  wave¬ 
lengths  using  a  dichroic  mirror.  Each  beam  is  then  focused 
onto  a  rotating  diffraction  grating  which  splits  each  beam  into 
several  pairs.  The  two  first-order  beams  for  each  wavelength 
are  then  recombined  onto  (lie  optical  axis  using  a  dichroic 
mirror,  collimated  and  focused  at  a  point  to  form  the  two- 
component  probe  volume.  In  the  present  study,  the  transmit¬ 
ting  optics  utilized  a  500  mm  focal  length  lens.  The  receiving 
optics  were  located  30  degrees  off  axis  in  the  forward-scatter 
direction.  Light  was  collected  using  a  500  mm  focal-length 
lens  and  then  focused  onto  a  100  pm  by  1  mm  long  slit.  The 
collected  light  is  then  split  and  picked  up  by  four  photo¬ 
detectors.  Three  are  arranged  to  look  at  the  signals  from  the 
514.5  beams  and  one  receives  light  from  the  488  nm  beams. 
Each  of  the  three  photodetectors  for  the  green  beams  are 
imaged  at  a  different  area  of  the  collection  lens  and  the  phase 
difference  between  the  signals  is  used  for  the  size  determina¬ 
tion.  Details  of  the  instalment  can  be  found  in  Ref.  7. 

In  the  present  study,  velocities  of  both  the  liquid  and  gaseous 
phases  were  measured.  This  was  accomplished  by  seeding  the 
gas  phase  with  nominal  I  pm  size  aluminuin-oxide  particles. 
The  coflow,  air-assist  flow,  and  the  ambient  surroundings 
were  all  seeded  to  minimize  biasing.  Phase  discrimination  is 
inherent  in  the  instrumentation  with  the  ability  to  size  each 
measured  particle.  At  each  spatial  location,  two  measurements 
were  taken  in  older  to  accurately  measure  the  velocity  of  each 
phase.  A  threshold  voltage  for  die  photodetectors  at  the  speci¬ 
fied  laser  power  was  determined  experimentally,  below  which 
signals  from  the  aluminum-oxide  particles  were  not  detected. 
For  the  droplet  measurements,  the  photodetcclor  voltage  was 
kept  below  this  threshold  value  m  order  to  eliminate  interfer¬ 
ence  from  the  aluminum-oxide  particles.  Total  laser  power  for 
all  wavelengths  was  fixed  at  1.5  W  for  ail  the  measurements. 
Particles  with  diameters  less  than  2.4  pm  were  used  to  repre¬ 
sent  the  gas  phase  velocity.  Two  complete  traverses  were 
taken  in  order  to  measure  all  three  components  of  velocity  and 
provide  a  check  on  flow  symmetry.  Each  traverse  measured 
axial  velocity  and  eidicr  radial  or  angular  velocity.  Generally. 
64  000  measurement  attempts  were  made  at  each  measurement 
location.  The  percentage  of  measurements  actually  validated 
depended  on  the  number  density  and  velocities  of  drops  at 
each  location  and  ranged  from  about  65  to  90  percent. 
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B.2  RESULTS  AND  DISCUSSION 

In  the  present  study,  results  are  presented  for  a  single  axial 
location  at  5  mnt  downstream  of  the  nozzle.  Gas  phase  results 
for  mean  velocities  are  presented  in  Figs.  3(a)  to  (c)  for 
isothermal,  single-phase  flow  without  droplets  and  two-phase 
flow  with  combustion.  Mean  gas  phase  axial  velocity, 
presented  in  Fig.  3(a).  presents  results  from  a  complete 
traverse  across  the  combustor  and  illustrates  the  symmetry  of 
the  flowfield.  The  combustor  exit  dimensions  are  illustrated  on 
the  x-axis  of  the  figure  for  reference.  A  small  recirculation 
zone  is  evident  near  the  center  of  the  nozzle.  At  this  axial 
location  relatively  close  to  the  nozzle,  velocity  gradients  are 
extremely  large  in  the  flow  from  the  air-assist  stream  contain¬ 
ing  the  droplets.  Effects  of  combustion  on  the  flowfield  are 
significant.  Both  the  maximun  and  minimum  mean  axial 
velocities  are  increased  for  the  combusting  case  compared  to 
the  isothermal  case  due  to  the  gas  expansion  associated  with 
the  heat  release.  For  the  case  with  combustion,  velocities 
increase  from  nearly  0  (o  35  m/s  and  then  decrease  to  -30  m/s 
within  a  radius  of  about  12  mm.  Mean  axial  velocities  in  the 
coflow  stream  are  not  affected  by  the  combustion  at  this  axial 
location. 

Figure  3b  presents  mean  radial  velocity  for  the  gns  phase.  As 
shown  in  Fig.  3(b),  effects  of  combustion  are  very  dramatic 
for  radial  velocity.  Maximum  radial  velocities  increased  from 
about  10  m/s  for  the  isothermal  case  to  about  40  itv's  for  the 
case  with  combustion  due  to  the  radial  expansion  of  the  gas. 
Again,  the  gas  from  the  coflow  stream  is  not  affected  by  the 
combustion  at  this  axial  location. 

Mean  gas  phase  angular  velocities  are  presented  in  Fig.  3(c). 
For  this  case,  reaction  and  the  presence  of  droplets  decreases 
the  maximum  angular  velocities  in  the  flowField.  Some  of  the 
decrease  in  angular  velocity  for  the  gas  phase  can  be  attri¬ 
buted  to  the  momentum  transferred  to  the  droplets  since  they 
do  not  initially  have  a  swirl  component. 

Fluctuating  gas  phase  velocities  are  presented  in  Ftgs.  4(a)  to 
(c)  for  both  (lie  single-phase,  isothermal  and  the  two-phase, 
combusting  cases.  All  fluctuating  velocities  presented  are  root- 
mean-squnred  (rms)  values.  Figure  4(a)  presents  radial  profiles 
of  gas  phase  fluctuating  axial  velocity.  Maximum  values  of 
fluctuating  axial  velocity  are  similar  for  both  the  combusting 
and  isothermal  cases.  The  case  with  combustion  does  show 
larger  values  of  axial  rms  velocity  at  radial  locations  between 
approximately  5  and  15  mm  from  the  center  of  the  nozzle. 
Axial  velocities  are  also  higher  at  these  locations  for  the 
combusting  case,  see  Fig.  3(a).  Fluctuating  radial  velocities, 
illustrated  in  Fig,  4(b).  show  dramatic  differences  between  (lie 
isothermal  and  combusting  cases.  The  maximum  velocity  loca¬ 
tions  have  shifted  radially  outward  corresponding  to  the  shift 
in  mean  radial  velocity,  see  Fig.  3(b).  The  maximum  fluctuat¬ 
ing  radial  velocity  has  also  increased  from  about  10  itv's  to 
15  m/s. 

Fluctuating  angular  gas  phase  velocities  are  presented  in 
Fig.  4(c).  Similar  to  the  results  shown  for  mean  angular  gas 
phase  velocities,  fluctuating  angular  velocities  generally 
decreased  with  combustion  and  the  presence  of  the  liquid 
phase  compared  to  the  single-phase,  isothermal  case.  A  small 
region  from  a  radius  of  about  7  to  15  mm  shows  increased 
values  of  fluctuating  angular  velocity  for  the  case  with  combustion. 


Mean  velocities  for  the  drops  are  presented  in  Figs.  5(n)  to  (c) 
for  llte  case  with  combustion.  In  the  experimental  study, 
velocities  were  measured  for  drop  sizes  ranging  from  4  to 
142  pm.  Results  are  presented  for  drop  sizes  of  15.  32  and 
52  pm.  Measured  gas  phase  velocities  are  also  presented  in 
the  figure.  Note  that  results  are  only  illustrated  from  -15  to 
+  15  mm  for  the  radial  direction  because  no  drops  were  present 
at  larger  radial  locations.  Figure  5(a)  presents  mean  drop  axial 
velocity  at  5  min  downstream.  Similar  to  the  results  previously 
shown  for  the  gas  pltase.  the  flowfield  is  very  symmetric. 
Axial  velocity  is  correlated  with  drop  size  in  all  regions.  In 
the  main  region  of  the  spray,  at  a  radius  of  about  7  mm,  the 
maximum  velocity  of  the  gas  phase  was  about  38  m  s.  and 
about  28  m/s  for  Ihe  32  pm  drops.  Even  the  maximum  veloc¬ 
ity  of  the  15  pm  drops  lagged  the  gas  pltase  by  about  5  m/s. 
In  the  center  of  the  flowfield  is  a  small  recirculation  zone,  see 
Fig.  3(a).  There,  only  die  15  pm  drops  showed  negative  axinl 
velocities  while  larger  drops  had  posiltve  velocities. 

Mean  drop  radial  velocities  are  presented  in  Fig.  5(b).  Again, 
there  is  a  correlation  between  drop  size  and  velocity  in  the 
flowfield.  Maximum  mean  radial  velocities  are  slightly  higher 
than  maximum  axial  velocities  for  the  drops  due  to  the  heat 
release  and  radial  expansion  of  the  gas.  Mean  angular  veloci¬ 
ties  of  the  drops  are  presented  in  Fig.  5(c).  Angular  velocity 
is  not  as  symmetric  and  is  also  much  smaller  than  the  other 
two  components  of  velocity.  The  mean  drop  angular  velocity 
is  a  strong  function  of  the  drop  size  with  the  smaller  drops 
showing  the  least  velocity  difference  with  the  gas  phase. 

Fluctuating  droplet  axial,  radial,  and  angular  velocity  com¬ 
ponents  for  the  three  drop  sizes  and  gas  phase  are  presented 
in  Figs.  6(a)  to  (c).  respectively.  The  fluctuating  drop  veloci¬ 
ties  presented  are  root-mean-squared  (rms)  values.  Generally, 
the  smaller  drops  are  affected  more  by  the  gas  phase  turbu¬ 
lence  and  have  larger  fluctuating  velocities  than  the  larger 
drops.  Velocity  fluctuations  are  clearly  not  isotropic  since 
fluctuating  axial  and  radial  velocities  are  considerably  larger 
than  fluctuating  angular  velocities. 

In  addition  to  drop  mean  and  fluctuating  velocities,  the  liquid 
volume  flux  is  important  in  two-phase  flows.  Drop  numher- 
flux  measurements  are  presented  in  Fig.  7,  where  results  are 
illustrated  for  four  drop  size  groups.  As  shown  in  Fig.  7. 
number  flux  is  nearly  symmetric.  The  results  show  that 
smaller  droplets  have  much  larger  number  fluxes.  The  distri¬ 
bution  of  the  larger  droplets  is  still  very  important  since  much 
of  the  liquid  mass  is  contained  in  the  Inrger  droplets. 
Relatively  few  drops  are  found  in  the  center  region  of  the 
flowfield  due  to  the  45  degree  s wirier  that  is  used  in  the  air- 
assist  stream. 

C.  COMBUSTION  CFD  CODE 

The  objective  of  the  present  work  is  to  develop  a  numerical 
solution  procedure  which  can  efficiently  handle  the  coupling 
between  a  spray  model  and  a  well-developed  strongly  implicit 
flow  solution  algorithm  (Refe,  13  and  14).  In  the  past,  spray 
models  have  been  coupled  with  different  flow  algorithms  and 
shown  some  promising  results  (Refs.  15  to  17).  However, 
most  of  Ihe  spray  models  were  coupled  with  a  flow  solver 
employing  a  segregated  approach,  such  as  a  TEACH-type 
code  (Refs.  15  and  16).  which  has  been  used  very  extensively 
in  the  industries  for  the  past  two  decades.  Although  simple 
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and  easy  to  implement,  the  TEACH-type  code  usually  suffers 
poor  convergence  due  to  the  explicit  (or  seini-implicit) 
treatments  of  the  chemical  source  terms  and  the  sequential 
solution  approach.  Recent  development  of  CFD  techniques  and 
the  advent  of  computer  technology  have  allowed  us  to  explore 
more  ambitious  schemes  to  solve  reacting  flow  problems. 
Strongly  coupled  and  implicit  numerical  schemes,  although 
requiring  much  more  computer  storage  and  complexity  of  the 
algorithm,  have  been  very  popular  for  nonreacting  compressi¬ 
ble  (Refs.  18  and  19)  and  incompressible  (Refs.  20  and  21) 
flow  computations. 

Shuen  and  Yoon  (Ref.  22)  developed  a  coupled  scheme  for 
high  speed  reacting  flows.  RPLUS.  which  has  been  used  and 
studied  quite  extensively  in  recent  years.  However,  like  other 
compressible  flow  solution  algorithms.  RPLUS  may  not  be 
suitable  for  low  s|>eed  flow  computations.  There  are  two  well- 
recognized  reasons  (Refs.  23  to  25)  for  the  convergence 
difficulties  related  to  compressible  flow  codes.  Fust,  the 
system's  eigenvalues  become  stiff  at  low  flow  velocities. 
Second,  the  pressure  term  in  the  momentum  equation  becomes 
singular  as  die  Mach  number  approaches  zero,  which  yields  a 
large  roundoff  error  and  smears  the  pressure  variation  field. 
This  will  not  only  result  in  slow  convergence  but  often  will 
produce  inaccurate  solutions.  To  circumvent  the  above  diffi¬ 
culties.  Shuen  et  al.  (Ref.  13)  developed  a  coupled  numerical 
algorithm  for  chemical  nonequilibrium  viscous  flows. 
ALLSPD.  which  utilizes  die  decomposition  of  die  pressure 
variable  into  a  constant  reference  pressure  and  a  gauge 
pressure  to  reduce  the  roundoff  errors  and  adds  a  precondition¬ 
ing  time  derivative  term  to  rescale  dte  system  eigenvalues. 
The  results  of  diese  treatments  show  diat  the  convergence 
properties  are  almost  independent  of  the  flow  Mach  number. 

Numerous  spray  models  have  been  proposed  and  investigated 
for  different  spray  combustion  problems  (Refs.  26  and  27)  in 
the  past  decade.  Recent  spray  models  differ  in  specific  details, 
but  generally  may  be  divided  into  two  categories:  locally 
homogeneous  flow  (LHF)  models  and  separated  flow  (SF) 
models.  LHF  models  represent  the  simplest  treatment  of  a 
multiphase  flow  and  have  been  widely  used  to  analyze  sprays 
(Ref.  28).  The  key  assumption  of  the  LHF  model  is  that  inter¬ 
phase  transport  rates  are  fast  in  comparison  to  the  rate  of 
development  of  the  flow.  This  implies  that  all  phases  have 
identical  properties  at  each  point  in  the  flow.  Clearly.  LHF 
models  are  only  formally  correct  for  flows  containing  infi¬ 
nitely  small  droplets. 

Numerous  SF  models  have  been  proposed  to  consider  inter¬ 
phase  transport  phenomena  (Ref.  29).  Among  them,  the  dis¬ 
crete  droplet  approach  (Refs.  26  and  30)  has  been  adopted, 
since  it  reduces  numerical  diffusion  while  providing  a  conve¬ 
nient  framework  for  dealing  with  multiple  droplet  size  and 
complex  uiterphase  transport  phenomena.  Many  discrete  drop¬ 
let  models  neglect  the  effects  of  turbulence  on  interphase 
transport  (Refs.  26  and  30).  This  implies  that  droplets  follow 
deterministic  trajectories,  yielding  the  deterministic  separated 
flow  (DSF)  model.  Neglecting  the  effects  of  turbulence  on 
droplet  transport  is  appropriate  when  characteristic  droplet 
relaxation  times  are  large  in  comparison  to  characteristic  times 
of  turbulent  fluctuations.  Few  practical  sprays,  however, 
satisfy  the  condition.  Dukowicz  (Ref.  31)  and  Gosnuin  and 
lonnnides  (Ref.  32)  have  adopted  stochastic  methods  to  study 
droplet  dispersion  by  turbulence.  Faeth  and  coworkers 


(Ref.  33)  extended  the  analysis  of  Gosman  et  al.  to  include  the 
effects  of  turbulence  on  interphase  heat  and  mass  transport. 
Their  stochastic  separated  flow  (SSF)  model  (Ref.  33)  has 
been  evaluated  in  a  wide  variety  of  parabolic  flows  with  very 
encouraging  results. 

In  a  recent  paper  (Ref.  34).  the  ALLSPD  algorithm  was 
extended  to  include  a  SSF  spray  model  (Ref.  33).  a  recently 
developed  low  Reynolds  number  k-e  turbulence  model 
(Ref.  35)  and  a  multiblock  treatment  to  calculate  (he  gas 
turbine  combustion  flows,  where  the  liquid  spray  is  an  im¬ 
portant  ingredient  of  combustion.  Although  the  turbulence 
model  has  been  included  in  this  study  for  turbulent  combus¬ 
tion  flow  calculations,  the  mean  flow  quantities  are  still  used 
in  the  chemistry  calculations.  The  turbulence  closure  problem 
for  chemistry  is  more  complicated  and  computationally  inten¬ 
sive.  The  consideration  of  a  suitable  turbulent  combustion 
closure  model  is  left  for  our  next-phase  study  in  the  near 
future. 


In  this  paper,  in  addition  to  the  gas  turbine  spray  combustion 
flow,  a  single-phase  premixed  turbulent  combustion  flow  and 
a  nonreacting  turbulent  flow  are  also  included  to  demonstrate 
the  current  status  of  this  research.  In  the  following  sections,  a 
brief  mathematical  formulation  of  the  governing  equations  for 
both  gas  and  liquid-phases  is  described.  The  numerical  method 
and  discretization  procedure  are  given  next  and  finally  some 
sample  numerical  results  are  presented. 


C.l  GOVERNING  EQUATIONS 

C.la  Gas-Phase  Equations 
C.lafl)  Navier-Stokes  Formulation 

The  two-dimensional,  unsteady,  compressible,  density- 
weighted  time-averaged  Navier-Stokes  equations  and  species 
transport  equations  for  a  chemically  reacting  gas  of  N  species 
written  in  generalized  nonorthogonal  coordinates  can  be 
expressed  as 
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where  the  vectors  Q,  E ,  F,  Ev,  Fv,  He,  and  H t  are  defined  as 
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In  the  above  expressions,  r.  5.  and  q  are  the  time  and  spatial 
coordinates  in  the  generalized  coordinates  and  5t  and  q,  are 
the  grid  speed  terms.  The  5*-  5,.  Hr  and  qv  are  the  metric 
terms  and  the  J  is  the  transformation  Jacobian.  The  power.  6. 
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is  an  index  for  Iwo  types  of  governing  equations  with  5  =  0  and 
for  two-dimensional  and  b  =  I  for  axisymmetric  cases  (with 
x  being  the  axial  and  v  the  radial  coordinates,  respectively). 

The  vectors  Q.  E.  K.  E^  Fv,  Hc  and  H(  in  the  above  defini¬ 
tions  are 


Q  =  (p,  p/i.  pv,  p £.  pK,  pc,  pK|. 


F,  -  (p/i,  p n1  ♦  /»,  p/i v,  (p£  +  /i)n.  puK. 

“  P"t.  P"Y ! . P',KN-l)  - 

F  "  (pv.  pnv.  pv:  +  /».  (p£  +  />)v,  pvK. 

"  pvt.  prF, . PvFn.,)7-. 


'V"V 

^  .  4it  J  ,lup 

i  p  t*  3  '/’>  r  llt 

r-.  4/T  3 

2^,  "  —  prri'nr  j 

Y.r  ”P'hrhft  "  47rr/T,V/'  Ar 

0 
0 

£,  vV 


Ev  =(o.tia.  T„.  HT,,  ♦  VT„  +  </„.  rA11. 

"  Tvt  ‘III . fl.  N-l)  • 


Fv  -(o.  Tn.  T„.  //rn  +  VTvv  ♦  l/vf.  T,*. 

*  Tv,  </,./ . fly  N-l 


and  the  source  term  vectors  Hc,  and  H,  are 
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where  p.//.  //.  v.  k.  c.  and  Y,  represent  the  density,  pressure. 
Cartesian  velocity  components,  turbulent  kinetic  energy, 
dissipation  rate  of  turbulent  kinetic  energy  and  species  mass 
fraction,  respectively:  £  =  e  +  *i(u2  *  v2)  is  the  total 
internal  energy  with  e  being  the  thermodynamic  internal 
energy;  and  5,  is  the  rate  of  change  of  species  i  due  to 
chemical  reactions.  The  normal  and  shear  stresses,  energy, 
species,  and  turbulent  diffusion  fluxes  are  given  by 
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where  7.  p^  pr  p,.  and  Jt,  are  (lie  temperature,  effective  vis¬ 
cosity.  turbulent  viscosity,  molecular  viscosity  and  the  effec¬ 
tive  thermal  conductivity,  respectively; 

"  <1  ~  */D,j 

is  the  effective  binary  diffusivity  of  species  i  in  the  gas 
mixture.  A,  the  molar  fraction  of  species  i,  and  Dt)  the  binary 
mass  diffusivity  between  species  i  and  j.  The  quantities  related 
to  the  source  term  in  the  turbulent  equations  are  given  as 
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is  the  specific  heat  of  the  gas  mixture  and  Pr,  is  the  turbulent 
Prnndtl  number. 

The  vector  H,  in  Eq.  (I)  represents  the  source  term  that 
accounts  for  the  interactions  between  the  gas  and  liquid 
phases.  In  vector  H,.  np  is  the  number  of  droplets  in  the  />th 
characteristic  group  of  droplets;  in.  is  the  evaporation  rate 
of  each  particle  group;  pp  is  the  liquid  density;  rp  is  the 
droplet  radius:  Up  and  vp  are  the  particle  velocities  at  the  /Hh 
group;  hfr  is  the  enthalpy  of  fuel  vapor  at  the  droplet  surface 
and  It  AT  is  (he  convective  heat  transfer  between  two  phases. 
Detailed  discussions  about  the  liquid  phase  equations  will  be 
described  later.  The  temperature  and  pressure  are  calculated 
iteratively  fioin  the  following  equations 
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l-l  p  w 
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where  R„  and  Tlr,  are  the  universal  gas  constant  and  reference 
temperature  for  thermodynamic  properties,  and  IV,.  CV,.  /i,,/ifl” 
are  the  molecular  weight,  constant  pressure  specific  heat, 
thermodynamic  enthalpy,  and  heat  of  formation  of  species  i. 
respectively. 

In  reacting  flow  calculations,  the  evaluation  of  thermophysical 
properties  is  of  vital  importance.  In  this  paper,  the  values  of 
C,„.  and  u„  for  each  species  are  determined  by  fourth- 
order  polynomials  of  temperature,  as  described  in  Shuen 
(Ref.  36).  The  specific  heat  of  the  gas  mixture  is  obtained  by 
mass  concentration  weighting  of  individual  species.  The 
thermal  conductivity  and  viscosity  of  the  mixture,  however, 
are  calculated  using  Wilke's  mixing  rule  (Ref.  37).  The  binary 
mass  diffusivity  D„  between  species  /  and  j  is  obtained  using 
the  Chapman-Enskog  theory  (Ref.  37). 

C.Ir(2)  All-Mach-Number  Formulation 

As  noted  earlier  in  the  introduction  section,  the  two  main 
difficulties  that  render  the  compressible  flow  algorithms 
ineffective  at  low  Mach  numbers  are  the  roundoff  error  caused 
by  the  singular  pressure  gradient  term  in  the  momentum  equa¬ 
tions  (the  pressure  term  is  of  order  I/M2  while  the  convective 
term  is  of  order  unity  in  the  nondimensional  momentum  equa¬ 
tions)  and  the  stiffness  caused  by  the  wide  disparities  in 
eigenvalues.  To  circumvent  these  two  problems  regarding  the 
low  Mach  number  calculations,  following  the  approach  by 
Merkle  and  Choi  (Ref.  23),  Shuen  et  al.  (Ref.  13)  added  a 
time  preconditioning  term  to  rescale  die  system  eigenvalues 
and  decompose  the  pressure  variable  into  a  constant  reference 
pressure  part  and  a  gauge  pressure  part.  This  all-Mach-number 
formulation  has  been  extended  to  include  the  turbulent  and 
spra)  equations.  The  resulting  Navier- Stokes  equations  in  a 
conservative  form  are 
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where  v„  in  the  expression  of  R,  is  the  normal  distance  away  whefe  the  priI,1IllvC  variable  veclorQand  the  preconditioned 

from  the  wall.  kt  is  the  molecular  thermal  conductivity,  Cpw  matrix  f  are  given  as 
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where  T*  is  the  pseudotime.  P  is  a  parameter  for  reseating  the 
eigenvalues  of  the  new  system  of  equations: 
and/i  =  t  *  p/pis  the  specific  enthalpy  of  the  gas  mixture. 
The  definition  of  vectors  in  Eq.  (3)  is  identical  to  those  in 
Eq.  (1)  except  that  the  pressure  terms  in  the  momentum 
equations  are  replaced  by  gauge  pressure.  />..  The  derivation 
of  the  all-Mach-number  formulation  can  be  found  in  detail  in 
Shuen  et  al  (Refs.  13  and  14). 


C.lb.  Liquid-Phase  Equations 
C.Ib(l)  Droplet  Motiou  Equations 

The  liquid  phase  is  treated  by  solving  Lagrangian  equations  of 
motion  and  transport  for  the  life  histories  of  a  statistically 
significant  sample  of  individual  droplets.  This  involves 
dividing  the  droplets  into  n  groups  (defined  by  position, 
velocity,  temperature  and  diameter)  at  the  fuel  nozzle  exit  and 
then  computing  their  subsequent  trajectories  in  the  flow.  The 
spray  model  used  in  this  study  is  based  on  a  dilute  spray 
assumption  which  is  valid  in  the  regions  of  spray  where  (lie 
droplet  loading  is  low  (Refs.  17  and  33).  The  liquid  fuel  is 
assumed  to  enter  the  combustor  as  a  fully  atomized  spray 
comprised  of  spherical  droplets.  The  present  model  does  not 
account  for  the  effects  due  to  droplet  breakup  and  coalescence 
processes  which  might  be  significant  in  a  dense  spray  situa¬ 
tion.  The  Lagrangian  equations  governing  the  droplet  motion 
are 
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The  subscript  x  represents  the  gas-phase  quantities  and  /; 
represents  the  liquid-phase  (or  'particle")  quantities.  Equations 
(4)  and  (5)  are  used  to  calculate  the  new  droplet  positions  and 
Eqs.  (6)  and  (7)  are  used  to  update  the  new  droplet  velocities 
at  the  new  droplet  locations.  A  second-order  Runge-Kutta 
scheme  was  used  to  integrate  Eqs.  (4)  to  (7). 

C.lb(2)  Droplet  Mass  and  Ueat  Transfer  Equations 

As  described  by  Faeth  (Ref.  26).  the  following  coiTelalions 
were  used  to  approximate  the  mass  and  heat  transfer  coeffi¬ 
cients  for  a  single  isolated  spherical  droplet: 
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where  h  is  the  heat  transfer  coefficient,  </  is  the  droplet 
diameter.  Dj  is  the  fuel  mass  diffusivity.  L  is  the  thermal 
conductivity  of  fuel  vapor  and  /hi  is  the  fuel  mass  evapora¬ 
tion  rate  per  unit  area.  The  Nr  araJ  are  defined  as 
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where  Sc  and  Le  are  the  Schmidt  and  Lewis  numbers,  respec¬ 
tively.  The  Spalding  number.  8.  in  Eqs.  (8)  and  (9)  is  defined 
as 
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where  the  particle  Reynolds  number.  Ri’n.  and  the  drag 
coefficient.  Cp  are  defined  as 


where  Y,<r  is  the  fuel  vapor  mass  fraction  at  the  surface  of  the 
droplet  nndF,f  is  the  mean  fuel  mass  fraction  of  the  ambient 
eas.  In  Hie  present  study.  Yltr,  is  obtained  from  the  following 
equation 
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where  W0  is  the  molecular  weight  of  gas  excluding  fuel  vapor. 
I Vf  is  the  molecular  weight  of  fuel  and  Xjgp  is  the  mole 
fraction  of  fuel.  The  X^p  is  obtained  from  the  assumption  of 
Raoult's  law.  Based  on  this  assumption,  the  mole  fraction  at 
the  droplet  surface  is  equal  to  the  ratio  of  the  partial  pressure 
of  fuel  vapor  to  the  total  pressure.  For  the  present  spray 
calculation,  the  partial  pressure  of  fuel  vapor  was  computed 
based  on  the  following  empirical  correlation  (Ref,  37): 
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where  x  -  I  -  7VT(.,  Pv,  *  -7.28936.  Pv2  -  1.53679. 

*  -3.08367.  Pv4  »  -1.02456.  7C  and  pc  are  tlte  critical 
temperature  and  critical  pressure  of  the  fuel  vapor, 
respectively. 


where  the  value  of  a  is  between  0  and  1  with  a  -  0  referring 
to  the  vortex  center  and  a*l  referring  to  the  droplet  surface. 
The  initial  and  boundary  conditions  for  Eq.  (13)  are 
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is  obtained  from  the  energy  balance  at  the 
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droplet  surface  by  the  following  equation 


C.lb(3)  Droplet  Internal  Temperature  Equations 


As  a  single  droplet  enters  a  hot  environment,  the  immediate 
small  portion  of  the  droplet  near  the  surface  will  be  heating  up 
quickly  while  the  center  core  of  the  droplet  remains  "cold". 
The  heat  will  be  conducted  and  convected  to  the  entire  interior 
as  the  droplet  penetrates  further  into  the  hot  ambient  gas. 
Eventually,  the  temperature  within  the  droplet  will  become 
nearly  uniform  before  the  end  of  its  lifetime.  To  solve  this 
transient  phenomena  within  the  droplet  is  not  trivial.  In  the 
past,  certain  approximations  (Ref.  38)  are  usually  made  to 
alleviate  this  computational  burden  while  obtaining  reasonably 
good  results.  Among  them,  the  simplest  one  is  die  uniform 
temperature  model.  This  model  assumes  that  the  thermal  con¬ 
ductivity  of  the  fuel  is  infinite.  Of  course,  this  is  not  valid  at 
the  beginning  of  the  heating-up  process  of  the  droplet. 
Another  model  considers  the  heat  diffusion  inside  the  droplet. 
Tins  yields  a  conduction  model.  The  temperature  distribution 
within  the  droplet  is  obtained  by  solving  the  one-dimensional 
heat  conduction  equation  subject  to  the  convective  boundary 
conditions  at  the  droplet  surface.  The  conduction  model 
completely  neglects  the  convective  phenomena  within  the 
droplet  which  might  occur  due  to  significant  shear  forces  at 
the  surface  induced  by  high  slip  velocities.  Tong  and 
Sirignano  (Ref.  39)  developed  a  vortex  model  for  the  internal 
temperature  of  a  single  droplet  which  accounts  for  the 
convective  effect  of  the  Hill's  vortex  formation  inside  the 
droplet.  In  this  study.  Tong  and  Sirignano’s  vortex  model  is 
applied  to  obtain  the  internal  temperature  distribution  of  the 
droplet.  The  equation  governing  the  internal  temperature 
distribution  based  on  this  model  is 
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where  and  It  are  calculated  from  Eqs.  (8)  and  (9).  respec¬ 
tively.  Iifg  is  the  latent  heat  of  the  fuel  and  AT  •  fg  -  Tpr . 
where  T/It  is  droplet  surface  temperature  and7^ is  the  mean  gas 
temperature  evaluated  in  the  following  way 
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An  implicit  scheme  was  used  to  solve  Eq.  (13)  subject  to  the 
initial  and  boundary  conditions.  Second-order  central  differ¬ 
ences  were  used  for  the  spatial  differential  terms  and  a  first- 
order  difference  for  the  time  term.  These  treatments  rendered 
a  scalar  tridiagonal  algebraic  system  and  was  solved  by  the 
Thomas  algorithm  (Ref.  40). 


C.2  NUMERICAL  METHODS 
C.2a  Discretized  Equations 

Equation  (3)  is  the  final  gas-phase  governing  equation  to  be 
solved  numerically.  To  obtain  time-accurate  solutions  for 
time-evolving  problems,  a  dual  time-stepping  integration 
method  can  be  applied  to  Eq.  (3).  The  solution  converged  in 
pseudotime  corresponds  to  a  time-accurate  solution  in  physical 
time  (Ref.  13).  However,  for  the  present  study,  since  only  the 
steady  state  solution  is  of  our  interest,  the  physical  time  term 
in  Eq.  (.3)  can  be  dropped  and  the  solution  can  be  marched 
completely  in  pseudotime  to  obtain  the  final  steady  state 
solution.  One  advantage  of  marching  to  the  steady  state 
solution  in  pseudotime  is  that  the  convergence  of  the  marching 
(iterative)  process  is  determined  by  the  eigenvalue  characteris¬ 
tics  on  the  pseudotime  space  and  not  by  the  original  stiff 
eigenvalues.  The  analysis  of  the  eigensystem  has  been  per¬ 
formed  mi  our  previous  study  (Refs.  13  and  14)  and  therefore 
will  not  be  repented  here.  It  should  be  noted  dial  the  inclusion 
of  the  turbulent  k-e  equations  in  this  study  does  not  affect 
the  system  eigenvalues  at  all  and.  therefore,  the  properties  of 
the  all-Mach-number  formulation  analysis  in  previous  papers 
(Refs.  13  and  14)  are  still  valid. 
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After  linearization  and  applying  a  first-order  time  differencing. 
Eq.  (3)  can  be  expressed  in  the  following  form: 
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Since  a  multiblock  treatment  is  applied  to  the  present 
numerical  algorithm,  the  interfaces  between  blocks  (or  zones) 
become  a  special  kind  of  boundary.  At  the  interface,  the 
governing  equations  are  still  solved  without  any  distinction 
from  the  rest  of  the  interior  points.  The  flow  variables  have 
been  carefully  arranged  in  such  a  way  that  the  information 
from  the  neighboring  biock(s)  is  automatically  brought  into 
the  calculation  at  the  interface  of  the  cunem  block.  The 
multiblock  treatment  of  the  present  study  not  only  provides 
flexibility  for  complex  geometry  calculations  but  also  reduces 
the  size  of  the  storage  array  for  the  MSIP  coefficients  due  to 
the  smaller  sizes  of  the  subdomams. 
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where  //  denotes  the  previous  iteration  level.  D  is  the  Jacobian 
for  chemical  and  turbulent  source  terms.  A  and  B  arc  die 
inviscid  term  Jacobians  and  Rcc  and  Ron  are  the  viscous  term 
lacobians.  The  expressions  for  these  Jacobians  can  he  found 
in  Shuen  et  al.  (Refs.  13  and  14)  except  for  the  turbulent  part. 
Central  differences  were  used  to  discretize  the  spatial  deriva¬ 
tive  terms  in  Eqs.  (15)  and  ( 16)  for  bodi  explicit  and  implicit 
operators.  The  resulting  coupled  algebraic  equations  are  solved 
using  a  modified  strongly  implicit  procedure  (MSIP)  proposed 
by  Schneider  and  Zedan  (Ref.  41 ).  which  is  completely  veclor- 
izable  along  /  +  2/  diagonal  direction. 

C.2b  Boundary  Conditions 

The  boundary  conditions  for  gas  phase  equations  arc  described 
as  follows.  Here,  only  the  subsonic  flow  boundary  conditions 
are  considered  in  the  present  paper.  At  the  inlet,  all  quantities 
are  specified  except  pressure  which  is  obtained  through  extra¬ 
polation  from  the  pressure  at  interior  points.  At  the  exit,  the 
governing  equations  are  solved  at  the  exit  station  by  applying 
backward  differences  for  the  slreamwise  derivative  terms 
(central  differences  are  still  used  for  cross-stream  derivative 
terms).  However,  flic  streamwise  pressure  derivative  terms  are 
i ■entrally-tUfferettced.  This  treatment  requires  the  pressure 
information  at  the  station  one  step  downstream  (outside  the 
computational  domain)  where  a  constant  pressure  condition  is 
enforced.  At  the  symmetry  line,  the  governing  equations  are 
solved  via  the  use  of  the  symmetry  conditions  for  two- 
dimensional  flows.  For  the  axisyininctric  case,  a  singularity 
exists  at  this  line  and  a  simple  one-sided  difference  is  used  to 
implement  the  symmetry  conditions.  At  the  solid  wall,  no  slip 
conditions  are  used  for  the  velocities  and  an  adiabatic  wall  is 
assumed.  Normal  derivatives  for  species.  Yf.  and  (j>f  +  2/3  pM 
are  set  to  zero.  For  the  turbulent  quantities,  the  low  Reynolds 
number  turbulence  model  (Ref.  35)  used  in  this  study  requires 
the  specification  of  the  k  and  c  at  the  walls  as  follows: 

k  =  0250m" 

4 

ft? 

c  =  0.251  _L 

v 

where  mt  is  the  friction  velocity  at  the  wall. 


C.2c  Spray  Source  Terms 

In  Eq.  (3).  the  liquid-phase  interaction  with  the  Navier-Stokes 
equations  is  modeled  as  a  source  term  which  accounts  for 
mass,  momentum  and  energy  exchanges  between  two  phases. 
In  order  to  compute  this  source  term,  the  liquid-phase  govern¬ 
ing  equations  described  in  the  previous  section  are  integrated 
in  time  once  initial  conditions  have  been  specified.  As  a  liquid 
droplet  begins  its  journey  from  the  injection  nozzle,  its 
position,  velocity,  temperature  and  size,  in  general,  will  be 
changing  according  to  (lie  interaction  between  the  gas-phase 
solution  and  the  spray  properties.  In  contrast  to  the  Lulerian 
approach  for  the  gas-phase  equations,  the  Lagntngian  treat¬ 
ment  for  the  liquid-phase  equations  requires  interpolation  of 
the  flow  quantities  from  the  Eulerian  grid  to  the  particle 
positions  and  redistribution  of  the  spray  source  terms  from  the 
particle  positions  to  the  Eulerian  grid  for  the  gas-phase.  Since 
the  time  step  for  die  spray  equations  is  usually  much  smaller 
than  that  for  the  flow  equations  (especially  for  the  present 
MSIP  scheme  with  which,  in  general,  a  large  CFL  number  can 
be  obtained)  and  a  large  number  of  particle  groups  are 
required  to  statistically  represent  the  spray  behavior,  it  is  very 
important  to  keep  the  numerical  efficiency  of  the  present 
implicit  scheme  from  being  severely  degraded  down  by  the 
coupling  (or  inter-action)  between  the  gas  and  liquid-phases. 
A  strategy  to  overcome  this  problem,  at  least  for  steady  state 
flows,  has  been  successfully  applied  to  the  present  study,  It 
will  be  described  in  this  section. 

C.2c(l)  Stochastic  Process 

As  mentioned  in  the  introduction  section,  there  are  two  types 
of  separated  flow  models  for  spray  computations.  One  is  the 
detenninislic  separated  flow  (DSF)  model  and  the  other  the 
stochastic  separated  flow  (SSF)  model.  For  turbulent  flow 
calculations,  the  DSF  model  completely  neglects  the  disper¬ 
sion  effect  due  to  the  turbulent  motion  and.  therefore,  mean 
gas  flow  quantities  are  used  to  evaluate  the  right-hand-side  of 
Eqs.  (6)  and  (7).  In  the  present  study,  the  SSF  model  is 
applied  to  account  for  the  dispersion  effects  on  spray  charac¬ 
teristics  In  this  model,  the  gas-phase  velocity  fluctua¬ 
tions./^  and  i  t, are  generated  by  randomly  sampling  a 
Gaussian  probability  density  distribution  having  a  stnndnid 
deviation  of  (2k/3)1/_.  The  instantaneous  velocities  are  then 
used  to  evaluate  Eqs.  (6)  and  (7j.  A  fixed  number  of  samplin¬ 
gs  is  conducted  for  each  group  of  particles.  The  final  spray 
quantities  are  obtained  by  averaging  the  results  (source  terms, 
trajectories,. .  . ,  etc.)  among  the  total  samples.  Details  of  the 
SSF  model  can  be  found  in  Reft.  26  and  33. 
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C.2c(2)  Determination  of  Spray  Time  Step 

For  spray  combustion  calculations,  there  are  several  lime 
scales  involved  in  the  flow  field  which  can  differ  by  several 
orders  of  magnitude.  The  chemical  reaction  time  scale  is 
usually  very  small  compared  with  Lite  rale  of  evolution  of  the 
gas  flow.  This  is  also  true  for  an  evaporating  spray.  To 
accurately  calculate  the  particle  trajectories,  size  and  tem¬ 
perature.  the  integration  time  step  has  to  be  small.  This  is 
especially  severe  as  the  droplet  becomes  smaller  and  smaller 
toward  the  end  of  its  lifetime.  For  an  unsteady  problem,  the 
time  step  for  the  entire  system  will  be  controlled  by  the 
smallest  time  step.  For  a  steady  state  calculation,  however,  the 
time  step  for  chemical  reaction  is  usually  not  a  problem  in  the 
present  formulation  due  to  the  implicit  treatment  of  the 
chemical  source  term.  Our  experience  indicates  that,  for 
single-phase  combustion  calculations,  the  same  CFL  number 
usually  can  be  used  for  both  nonreacting  and  reacting  calcu¬ 
lations  with  the  present  MS1P  method.  This  ensures  that  the 
convergence  properties  for  reacting  flow  calculations  are  not 
degraded  using  the  present  numerical  algorithm.  However,  the 
time  step  for  spray  equations  still  remains  small  and  has  to  be 
selected  (computed)  very  carefully  in  order  to  obtain  accurate 
spray  results  and  to  maintain  stability.  The  determination  of 
the  spray  time  step  will  be  discussed  here  and  how  the  spray 
interacts  with  the  gas  flow  is  described  in  the  next  section. 

The  spray  time  step  at  any  instant  of  time  along  its  trajectory 
is  determined  based  on  (lie  following  time  step  constraints:  (1) 
droplet  velocity  relaxation  time  (fr).  (2)  droplet  life  time  (fy). 
(3)  droplet  surface  temperature  constraint  time  (rt).  (4)  local 
grid  time  scale  ( / g)  and  (5)  turbulent  eddy-droplet  interaction 
time  (I/).  The  final  spray  time  step  (A/r;lr)  is  determined  by 
taking  the  minimum  of  the  above  five  time  steps  to  ensure  the 
accuracy  and  stability  of  (lie  spray  calculations.  A  factor 
between  0.1  and  0.5  is  further  used  to  multiply  the  selected 
time  step  in  the  current  spray  calculation.  These  time  steps  are 
described  as  follows. 


In  order  to  ensure  the  success  of  the  temperature  calculation 
using  Eq.  (13).  the  exact  solution  of  (he  infinite  conductivity 
model  (or  unifonn  temperature  model)  is  used  to  estimate  the 
time  step  for  the  present  vortex  model.  The  temperature 
equation  for  the  droplet  based  on  this  model  is 


Equation  (17)  has  an  exact  solution  (after  local  linearization) 
in  the  following  form 

&Tp  =  a(]  -e 

where  A Tp  is  the  droplet  temperature  change  within  the 
integration  time  step.  A,  A',  ami  B'  are  defined  as 
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Therefore,  a  time  scale  can  be  obtained  based  on  the  above 
solution  if  a  desited  ATp  is  specified.  This  time  scale  is 
expressed  as 


In  our  spray  calculation.  AT^  “  3  K  is  specified  for  those 
particles  just  leaving  the  injector,  where  they  experience  a 
sudden  temperature  jump  and  A Tp  -  0.5  K  for  the  rest  of  the 
calculation  toward  the  end  of  their  lifetime. 


Droplet  velocity  relaxation  time  (I,.):  The  local  linearized 
droplet  equations  of  motion.  Eqs.  (6)  and  (7).  have  exact 
solutions  in  terms  of  the  local  slip  velocity  with  an  expo¬ 
nential  decay  fonn.  The  time  constant  for  the  exact  solution 
can  be  expressed  as 
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Droplet  fife  time  <t(>:  To  ensure  that  the  drop  size  remains 
positive  for  the  practical  computational  purpose,  the  droplet 
lifetime  at  any  instant  of  time  is  estimated  by  the  following 
equation 
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Droplet  surface  temperature  constraint  lime  (t,):  When  the 
governing  equation  for  the  droplet  internal  temperature 
distribution.  Eq.  (13).  is  solved,  the  temperature  solution  can 
become  completely  incorrect  due  to  the  use  of  an  inappro¬ 
priately  large  time  step.  This  is  particularly  important  for  fuel 
with  n  low  boiling  temperature  (close  to  room  temperature). 


Local  grid  time  scale  (1^):  A  particle  can  travel  across  several 
grids  and  may  experience  a  sudden  change  of  the  local  gas 
properties  if  the  time  step  is  loo  large.  This  not  only  causes 
inaccuracy  in  the  integration  but  also  increases  the  difficulties 
of  locating  the  particle  positions.  Therefore,  a  time  scale.  ig. 
is  computed  to  ensure  dial  the  particle  only  moves  less  titan 
one  local  cell  size  in  one  time  step. 

Turbulent  F.ddy-droplet  interaction  time  <t|):  According  to 
Shuen  et  al.  (Ref.  42).  a  particle  is  assumed  to  interact  with  an 
eddy  for  a  time  which  is  the  minimum  of  either  the  eddy  life¬ 
time  or  the  transit  time  required  for  the  particle  to  cross  the 
eddy.  These  times  are  estimated  by  assuming  that  the  charac¬ 
teristic  size  of  an  eddy  is  the  dissipation  length  scale  as 

,  -W 

K  “  K  /e 


and  the  eddy  lifetime  is  estimated  as 

ie  =  Le/(2  k/3)  1/2 . 

The  transit  time  of  a  particle  was  found  using  the  linearized 
equation  of  motion  for  a  particle  in  a  uniform  flow 


t,  -  -Tfn_l-  L,(t|  (T  -  iTp  |) 
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and  u  -  ii.  .  is  the  relative  velocity  nt  the  start  of  the 
interaction.  When  Lt  >  t\u  -  up  | ,  tlie  lineanzed  stopping 

distance  of  the  particle  is  smaller  than  the  characteristic  length 
scale  of  the  eddy  and  Eq.  (18)  has  no  solution.  In  this  CBse. 
the  eddy  hns  captured  the  particle  and  the  interaction  time  is 
the  eddy  lifetime.  Therefore 

‘t  m  <e-  •/  Lc  >  Tl«*  "  I 

f,  -  min (tt.  I,),  if  Le  <  t|/T  -  up  | 

C.2c(3)  Interaction  Between  Two  Phases 

For  the  gas-phase  equations.  Eq.  (3).  the  presence  of  the  spray 
appears  in  the  form  of  a  source  term.  H,.  This  source  term 
represents  the  interchange  of  the  mass,  momentum  and  energy 
between  two  phases.  As  the  particles  are  injected  into  the  flow 
domain,  their  subsequent  behavior  (positions,  velocities,  size, 
and  temperature)  is  affected  continuously  by  the  neighboring 
gas  properties  and  nee  versa.  Usually  a  large  number  of  spray 
particles  is  desirable  to  accurately  predict  the  spray  behavior. 
However,  this  requires  tremendous  computational  effort  for  the 
spray  calculation  alone.  To  minimize  the  computational  lime 
for  the  spray,  the  spray  source  term.  H(.  is  not  required  to  be 
updated  at  every  gas-phase  iteration  (pseudotime  time  march¬ 
ing).  Usually  the  spray  source  term  is  updated  every  10  to 
20  iterations  ui  our  spray  calculations.  When  the  spray  source 
term  is  updated,  each  group  of  particles  is  integrated  either  to 
the  end  of  its  lifetime  or  until  it  leaves  the  computational 
domain.  It  should  be  noted  that,  for  the  present  steady  state 
spray  combustion  calculation,  the  spray  tune  step  determined 
previously  is  independent  of  the  pseudotime  used  for  the  gas- 
phase  equations,  which  is  determined  mainly  from  the  obtain¬ 
able  maximum  CFL  number  accord-ing  to  the  local  system 
eigenvalues.  Therefore,  this  strategy  maintains  both  the 
efficiency  of  tlie  flow  solver  and  the  accuracy  of  the  spray 
Lagrangian  integration.  This  treatment  of  die  gas-liquid 
interaction  is  different  from  those  reported  by  Raju  and 
Sirignano  (Ref.  17)  where  time-accurate  solutions  were  their 
primary  concern.  In  die  present  computation,  it  is  assumed 
(lint  upon  impingement  with  die  walls,  (he  droplets  evaporate 
completely  and  assume  the  local  gas  flow  velocities.  The 
interpolation  of  the  gas-phase  properties  from  the  Eulerian 
grid  to  the  particle  location  and  the  redistribution  of  the  spray 
source  term  from  the  particle  location  to  the  Eulerian  grid  is 
applied  in  a  similar  way  as  discussed  in  Raju  and  Sirienano 
(Ref.  17). 

CJ  NUMERICAL  TEST  RESULTS 

In  this  section,  results  obtained  from  the  ALLSPD  algorithm 
with  and  without  spray  are  presented.  These  include  (1)  a 
nonreacting  turbulent  backward-facing  step  flow  to  demon¬ 
strate  die  validity  of  the  current  turbulence  model.  (2)  a 
single-phase  turbulent  reverse  jet  combustion  flow  to  assess 
the  present  combustion  treatment  and  (3)  a  spray  gas  turbine 


combustion  flow  to  qualitatively  demonstrate  the  spray 
calculation  for  a  complex  geometry  and  the  interaction 
between  die  two  phases.  The  spray  results  shown  here 
emphasize  the  numerical  aspect  radier  than  the  spray  physics. 
The  accuracy  validation  for  spray  calculations  will  be  consid¬ 
ered  in  future  calculations. 

C.3a  Backward-Facing  Step  Mow 

The  turbulent  backward-facing  step  flow  data  of  Kim  et  al. 
(Ref.  43)  for  a  two-dimensional  channel  with  an  inlet  to  step 
height  ratio  of  two  is  selected  here  to  test  the  validity  of  the 
k-c  turbulence  model.  A  136  *  100  grid,  clustered  near  the 
step  and  the  top  and  bottom  walls,  was  used.  No  chemical 
reactions  were  included  in  the  calculation. 

The  particle  traces  of  the  flow  is  shown  in  Fig.  8.  The 
experimental  reattachment  length  given  by  Kim  et  al.  is  7.1  II. 
where  H  is  the  step  height.  The  predicted  value  in  our  calcula¬ 
tion  is  about  6.1  H.  which  represents  a  14  percent  under- 
prediction.  Figure  9  shows  the  mean  velocity  profiles  al 
various  axial  locations.  The  agreement  is  quite  good  in  all 
locations  except  near  the  reattachment  point.  The  profiles  for 
the  turbulent  kinetic  energy  and  turbulent  shear  stress  are 
given  in  Fig.  10.  Good  agreement  is  observed  in  diesc  com¬ 
parisons.  Wall  pressure  coefficients  along  the  step-side  and 
opposite  walls  are  shown  in  Fig.  1 1 .  The  underprediction  of 
reattachment  length  is  also  evident  from  this  figure.  The  pre¬ 
dicted  pressure  recovery  downstream  of  the  reattachment  point 
is  in  excellent  agreement  with  measured  pressure  levels.  The 
convergence  history  for  tiiis  calculation  is  illustrated  in 
Fig.  12.  The  convergence  property  for  the  calculation  is 
satisfactory. 

C.3b  Reverse  Jet  Combustion  Flow 

The  flow  configuration  is  a  51  mm  I.D.  (inner  diam)  times 
457  mm  cylindrical  chamber  containing  a  reverse  jet  flame 
holder  which  issues  from  a  1.32  mm  I.D.  (6.35  mm  O.D.) 
tube.  The  jet  is  coincident  with  the  chamber  axis  and  located 
80  mm  upstream  from  the  chamber  exit.  Both  the  main  and  jet 
flows  are  stoichiometrically  preinixed  propane  and  air  at  a 
temperature  of  300  K.  with  a  mean  velocity  of  7.5  itvs  for  the 
main  stream  and  135  m/s  for  the  jet.  A  complete  description 
of  the  flow  system  is  available  in  McDannel  et  al.  (Ref.  44). 
A  1 17*01  grid  (half  domain)  clustered  near  the  jet  lube  was 
used.  Five  species  (Cyf^.  02.  N2,  C02.  and  H:0)  were 
considered  in  this  calculation  and  the  single-step  global 
reaction  chemistry  model  reported  in  Westbrook  and  Dryer 
(Ref.  45)  was  used  for  combustion.  Figures  13.  14.  and  15 
show  the  particle  traces,  velocity  vectors  (colored  by  tem¬ 
perature)  and  temperature  contours  of  the  reacting  flow, 
respectively.  The  experimentally  measured  temperature  con¬ 
tours  (directlv  taken  from  McDannel  et  al.  (Ref.  44))  are  also 
presented  in  Fig.  15  for  comparison.  The  flow  is  clearly  seen 
to  consist  of  two  distinct  regions  -  the  recirculation  zone  and 
the  wake.  The  incoming  flow  is  ignited  by  the  hot  combustion 
gas  in  the  recirculation  zone  and  further  combustion  takes 
place  in  tlie  wake.  For  the  test  conditions  considered  here  the 
reverse  jet  serves  as  a  very  effective  flame  holding  device. 
The  results  in  Fig.  15  indicate  that  the  predicted  temperatures 
are  higher  than  the  measured  values.  This  is  mainly  attributed 
to  the  over-simplified  chemistry  model  used  in  the  present 
calculation.  As  reported  by  McDannel  et  al.  (Ref.  44).  there 
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was  n  significant  amount  of  carbon  monoxide  (CO)  observed 
in  the  combustion  products,  which  lowered  the  combustion 
temperature.  Also  not  included  in  the  calculation  is  die  radia¬ 
tion  heat  loss  which  will  lower  the  flow  temperature  further. 
Ina  addition,  the  thermocouple's  used  in  the  measurements 
were  not  corrected  for  radiation  loss,  indicating  that  the  actual 
gas  temperatures  would  be  somewhat  higher  then  measured. 

C.,3c  Gas  Turbine  Spray  Combustion  Flow 

A  simplified  model  of  General  Electric's  EEE  (Energy 
Efficiency  Engine)  combustor  is  considered.  Figure  16  shows 
the  clustered  8  I  *65  grid  for  this  calculation  and  the  general 
engine  grid  layout  for  (he  entire  engine  (Ref.  46).  It  is  an 
nnnular  combustor  with  a  dual  can  combustor  dome.  The 
cooling  air  through  the  internal  walls  (two  combustor  domes) 
was  omitted  for  simplicity.  Also,  the  swirling  effect  was  not 
included  in  the  calculation.  The  present  calculation  does  not 
completely  simulate  the  typical  gas  turbine  combustion 
characteristics,  in  which  swirling  and  cooling  are  the  two 
important  ingredients.  We  would  like  to  emphasize,  in  the 
present  study,  that  the  focus  is  to  demonstrate  the  effective¬ 
ness  of  the  interaction  between  the  two  phases  based  on  the 
present  spray  solution  procedure.  Further  studies  with  this 
algorithm  will  be  conducted  both  on  the  spray  accuracy 
evaluation  and  on  the  detailed  swirling  and  cooling  computa¬ 
tion  for  gas  turbine  combustor  configurations. 

For  the  present  spray  combustion  calculation,  the  spray 
injectors  were  located  close  to  the  inlet  of  the  combustor 
domes.  One  hundred  spray  groups  with  10  random  samples  for 
each  groups  were  used,  which  can  be  seen  in  Fig.  17.  The 
case  studied  here  has  a  flow  Reynolds  number.  Re  m  1 .05*  Kr. 
where  the  Reynolds  number  is  based  on  the  inlet  maximum 
velocity  and  combustor  inlet  height.  The  inlet  air  temperature 
is  900  K  and  the  pressure  is  I  atm.  The  liquid  n-pentane  fuel 
was  used  and  five  species  (C]Ht2.  N2,  CO 2,  and  H20 ) 

were  considered  in  this  calculation.  Again,  the  single-step 
global  reaction  chemistry  model  reported  in  Westbrook  and 
Dryer  (Ref.  45)  was  used  for  combustion.  The  fuel/air  ratio  is 
0.02  (total  fuel  injected/total  incoming  air  at  the  inlet,  includ¬ 
ing  bypass  air).  The  injection  velocity  of  the  liquid  fuel  at  the 
exit  of  the  fuel  nozzle  is  20  rtv's  and  the  temperature  is  290  K. 
The  liquid  fuel  was  assumed  to  be  fully  atomized  with  the 
initial  diameters  ranging  from  20  to  100  pm.  The  liquid  fuel 
was  injected  into  the  gas  flow  after  the  gas  flow  had  been 
iterated  to  reach  a  nearly  steady  state  solution.  Upon  the 
injection  of  die  fuel,  the  cool  fuel  was  suddenly  exposed  to  a 
hot  environment  nnd  the  interaction  between  the  two  phases 
took  place  in  terms  of  the  interchange  of  mass,  momentum 
and  energy.  An  ignition  source  was  placed  downstream  of  die 
injector  to  ignite  the  "burning'  of  the  fuel-air  mixture.  The 
ignitor  was  turned  on  right  after  the  initiation  of  the  spray  and 
was  turned  off  when  the  temperature  in  any  of  the  ignition 
computational  cells  reached  1600  K.  The  spray  source  term. 
H|,  in  Eq.  (3)  was  updated  (computed)  every  20  gas-phase 
iterations.  The  converged  results  are  presented  in  the  following 
figures. 

Figure  10  shows  the  particle  trajectories.  The  dots  in  the 
figure  are  not  scaled  to  indicate  the  relative  size  of  the  actual 
liquid  particles,  although  smaller  particles  evaporate  much 
faster  titan  bigger  ones.  The  computed  velocity  vectors  colored 


with  temperature  are  presented  in  Fig.  18.  where  a  converged 
solution  of  the  cold  flow  (nonspray/noncombustion)  with  (he 
same  flow  condition  is  also  included  for  comparison.  The 
combustion  zones  are  confined  near  the  inner  walls.  This  is 
mainly  due  to  the  exclusion  of  the  swirlers  at  the  inlet  of  the 
combustor  domes  in  the  present  calculation.  Without  swirl,  the 
degree  of  the  fuel-air  mixing  is  relatively  poor.  As  may  be 
noted,  there  is  no  flame  holding  device  for  the  present 
calculation.  The  recirculation  zones  near  the  inner  walls  are 
the  only  devices  to  hold  the  flame.  Figure  19  shows  the  Sauter 
mean  diameter.  D}2.  along  the  averaged  trajectory  for  both  the 
lower  and  upper  domes.  The  increase  of  the  032  right  after  (he 
injection  indicates  the  rapid  evaporation  of  tire  smaller 
particles.  After  reaching  a  peak,  both  mean  diameters  gradu¬ 
ally  decrease  as  more  nnd  more  particles  evaporate.  To  show 
the  convergence  properties  for  both  nonreacting  and  spray 
combustion  cases,  the  convergence  histories  for  both  cases  are 
shown  in  Fig.  20.  As  can  be  seen  in  this  figure,  the  spray  was 
initiated  after  the  gas-phase  solution  had  been  iterated  1000 
times.  A  sudden  disturbance  from  the  spray  injection  causes 
the  L2  norm  residual  to  jump  to  a  level  higher  than  the  liutial 
gas-phase  residual.  As  more  and  more  spray  particles  evapo¬ 
rate  and  undergo  combustion,  the  interaction  between  two 
phases  can  be  clearly  identified  in  this  convergence  pattern. 
Since  in  the  present  calculation,  the  spray  source  term  was 
updated  every  20  gas-phase  iterations,  a  small  residual  spike 
along  the  convergence  history  can  be  seen  very  clearly.  These 
small  spikes  persist  toward  the  end  of  the  present  computation. 

D.  CONCLUDING  REMARKS 

Much  work  is  required  to  develop  a  CFD  code  to  the  level 
necessary  for  a  designer  to  be  able  to  use  it  with  confidence. 
We  at  NASA  are  proceeding  with  a  commitment  to  do  this.  A 
three-dimensional  version  of  our  ALLSPD  will  be  developed 
shortly  and  extensive  validation  against  experimental  data  will 
be  performed.  At  the  same  time  we  are  exploring  the  use  of 
advanced  modeling  to  better  represent  the  flow  physics  and 
chemistry  of  turbulent  combustion.  For  example,  modeling  of 
dense  sprays  and  sprays  at  super-critical  conditions  will  be 
investigated.  Modeling  of  the  turbulence  and  chemistry  inter¬ 
action  will  be  explored,  using  techniques  such  as  PDF 
methods.  The  use  of  reduced  hydrocarbon  chemical  kinetics 
models  which  represent  the  actual  combustion  processes  aie 
also  necessary.  And  finally,  an  improved  model  to  represent 
potentially  high  levels  of  radiation  heat  transfer  is  also  needed. 

The  state  of  art  for  numerics  also  continues  to  advance  at  a 
rapid  pace.  Our  code  will  he  modified  in  the  future  as  oppor¬ 
tunities  arise  in  areas  such  as  unstructured  grids  and  massively 
parallel  computing.  This  ongoing  effort  to  produce  a  modern 
chemical  reacting  flow  CFD  code  holds  promise  to  provide  a 
powerful  design  tool  for  the  industry's  use  in  the  analyses  of 
the  next  generation  of  gas  turbine  engine  combustor  concepts. 
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Figure  1 . — Schematic  drawing  of  the  combustor. 
Dimensions  in  mm. 


Figure  2. — Optical  configuration  of  the  phase/Doppler 
particle  analyzer. 


Mean  angular  velocity,  m/s  Mean  radial  velocity,  m/s  Mean  axial  velocity,  m/s 


15 


o  Single  phase,  isothermal 
•  Combusting  spray 


Radius,  mm 

(c)  Mean  angular  velocity. 

Figure  3. — Gas  phase  mean  velocity  at  5  mm  down¬ 
stream. 
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Figure  4. — Gas  phase  fluctuating  velocity  at  5  mm  down¬ 
stream. 
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Figure  5. — Drop  mean  velocity  at  5  mm  downstream. 
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Figure  6. — Drop  fluctuating  velocity  at  5  mm  downstream. 
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Figure  8. — Particle  traces  tor  turbulent  backward-facing  step  flow. 
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Figure  9. — Numerical  predictions  of  the  mean  velocity  profiles  at  selected  axial  locations  tor  the  turbulent 
backward-facing  step  flow. 
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(b)  Opposite  wall. 


Figure  1 1 . — Comparison  of  the  wall  pressure  coefficients 
for  the  turbulent  backward-facing  step  flow. 


Figure  13. — Particle  traces  for  the  reverse  jet  combustion 
flow. 


Figure  1 4. — Velocity  vectors  (colored  by  temperature)  for 
the  reverse  jet  combustion  flow. 


Figure  12.— Convergence  history  for  the  turbulent  backward¬ 
facing  step  flow. 
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(a)  Experimental  data  (from  McDannel  et  al.). 


Figure  1 5. — Comparison  of  the  temperature  predictions 
with  the  experimental  data  for  the  reverse  jet  combustion 
flow. 
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Figure  16. — Center-plane  cut  for  GE's  EEE  gas  turbine 
engine  and  the  grid  for  the  2-D  combustor. 


Figure  17.— Liquid  fuel  particle  trajectories  for  the  2-D 
gas  turbine  spray  combustion  flow. 


Figure  1 8. — Velocity  vectors  (colored  by  temperature)  for  (a)  non-combustion  and  (b)  spray  combustion  cases  for  the  2-D 
gas  turbine  spray  combustion  flow. 
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Figure  19.— Sauter  mean  diameter  distribution  for  upper 
and  lower  combustor  domes  along  the  averaged  tra¬ 
jectory  for  the  2-D  gas  turbine  spray  combustion  flow 


Figure  20. — Convergence  history  for  both  non-combus¬ 
tion  and  spray  combustion  cases  for  the  2-D  gas  tur¬ 
bine  spray  combustion  flow. 
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